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ABSTRACT 

this  paper,  the  principles  of  the  synthetic 
analysis  and  treatment  techniques  concerning  vehicle 
re-entry  thermal  environment,  ablation,  heat*  transfer 
and  thermal  stresses  are  described.  The  thermochemical 
ablation  calculations  for  chemical  reactions  on  the 
ablating  surfaces  of  carbon-base  materials  are  simpli¬ 
fied  using  selective  calculation  results.  The  temper¬ 
ature  fields  are  calculated  using  the  three-level 
explicit  difference  schemes  for  both  constant  and  non¬ 
constant  space  steps  within  the  total  region.  The 
thermal  stresses  are  calculated,  based  on  the  incomplete 
generalized  potential  energy  variational  principles 
under  the  supposed  conditions  of  linear  elasticity, 
using  the  finite  element  method.  This  method  can  also 
be  used  to  calculate  the  effects  of  the  various  displace 
ment  boundary  conditions.  An  automatic  grid  generation 
technique  has  been  developed  during  the  computer  code 
editing,  which  helps  improve  the  shape  selection  cal¬ 
culations. 
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1 .  INTRODUCTION 


It  has  long  been  a  subject  of  interest  to  study  and  develop 
vehicle  heat  protection  nosetips  using  carbon-base  materials  [1,3]. 
The  most  appealing  advantages  of  the  carbon-base  nosetips  are 
slight  ablation  and  small  shape  changes.  However,  the  thermal 
stress  they  experience  under  a  shock  heating  environment  is  a 
problem  of  major  concern.  In  order  to  investigate  the  nosetip 
heat  protection,  it  is  necessary  to  develop  a  computer  code  for  a 
synthetic  analysis  algorithm  to  simultaneously  calculate  the  thermal 
environment,  ablation,  heat  transfer  and  thermal  stresses.  The 
present  paper  only  gives  examples  for  the  plug-type  nosetipe  anal¬ 
ysis,  though  the  algorithm  discussed  here  can  also  be  applied  for 
an  analysis  of  shell-type  nosetips. 

For  a  heat  protection  nosetip,  the  state  of  the  thermal  stress 
depends  not  only  on  the  mode  of  its  structure,  but  also  on  the 
heat  environment,  ablation,  the  history  of  its  shape  changes  and 
the  instantaneous  temperature  fields  in  the  body.  Moreover,  all 
these  factors  interfere  with  each  other.  Every  effort  has  been 
made  in  recent  studies  to  fully  couple  these  factors  in  order  to 
simulate  realistic  conditions  as  closely  as  possible. 

«• 

2.  ABLATION  AND  ABLATING  CONTOUR 


During  re-entry,  the  thermochemical  ablation  of  the  heat  protec¬ 
tion  layers  of  the  carbon-base  materials  occurs  mainly  in  an 
extremely  high  tmeperature  range  which  is  over  3000  °K.  At  such 
high  temperatures,  the  thermochemical  reactions  occurring  on  the 
graphite  surface  can  be  considered  as  in  a  state  of  equilibrium 
[4] ,  During  re-entry,  the  temperature  on  the  surface  of  the  heat 
protection  layers  increases  quickly.  At  the  beginning  of  the 
re-entry,  there  exists  an  influence  of  the  non-equilibrium  chemical 
reactions.  Since  the  quantity  of  the  ablative  mass  loss  during 
this  period  is  limited,  it  should  not  introduce  appreciable  errors 
when  calculating  the  ablating  processes  during  the  whole  re-entry 
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period  under  the  assumption  of  chemical  equilibrium. 


Detailed  analysis  of  the  thermochemical  ablation  of  the  carbon 
base  nosetips,  based  on  the  JANAF  thermodynamics  data  [6],  and  the 
screening  calculations  on  the  entire  ablation  surfaces  along  the 
trajectory  for  five  carbon  vapors  and  10  carbon  compounds,  C  (gas) , 

^2'  ^3'  ^4'  ^5'  ^2°'  ^3°2'  CN2 (C-N-N) ,  CN2(N-C-N), 

*"2^'  *"3^2 '  ^4^2'  thermodynamic  data  are  available,  are  given 

in  [5] .  Results  at  the  stagnation  points  are  shown  in  Figure  1 
( i ' 

where  f^  '  (i  =  1 , 2 , 3 , . . . , 15)  are  the  fractions  of  the  ablated 
carbon  in  the  total  rate  of  thermochemical  ablation  in  the  chemical 
reactions  to  produce  each  gas  or  compound.  Results  show  that  on 
the  ablation  carbon  surfaces,  the  chemical  reactions  producing  CO, 
C,,  C^N  and  CN  (corresponding  to  curves  f^^^,  f^^^,  and 

in  Figure  1)  contributes  more  to  the  carbon  ablation  rate  and/or 
to  the  thermal  effects  than  to  other  effects.  This  makes  the  abla¬ 
tion  calculations  much  simpler.  Here,  some  of  the  secondary 
chemical  reactions  can  be  neglected  during  the  calculations,  and 
it  is  sufficient  to  take  into  account  only  the  principal  chemical 
reactions  and  corresponding  chemical  compositions,  i.e., 


3C=f=^C, 

C-t-l-O.^CO 

C  +  ^N,^tN 

2C+-|-N,^C.N. 


(1'  ) 


Thermodynamic  data  show  that  within  the  normal  temperature  range 
at  the  ablation  surface,  the  equilibrium  constants  and  the  thermal 
effects  of  these  reactions  can  be  expressed  nicely  by  the  fitting 
formulas  shown  in  Toble  1  and  the  approximate  mean  values.  This 
makes  the  ablation  calculations  very  convenient.  The  equilibrium 
constant  of  the  carbon-oxygen  reaction  can  be  replaced  by  zero 


of  the  oxygen  remainder.  Based  upon  its  equilibrium  constant,  the 
actual  quantity  of  the  oxygen  remainder  in  this  reaction  is  small 
enough  to  be  negligible. 

The  thermochemical  ablation  rate  of  mass  loss  on  the  wall 
surfaces  can  be  expressed  as 


rh^  =  - 

(1 ) 

where  the  dimensionless  rate  of  ablation,  B,can  be  found  through 
simultaneously  solving  the  gas  diffusivity  equations  in  the  boundar 
layer  with  the  equilibrium  constants  of  the  chemical  reactions. 

The  heat  flux  which  gets  into  the  burning  wall  is 

(2 ) 

9o..=  09o  +  gr,-£uri. 

TABLE  1.  The  approximate  equations  of  the  equilibrium  constants 
and  thermal  effects 
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where  the  second  and  third  terms  on  the  right  of  the  equation  are 
the  radiation  heat  of  the  air  to  the  wall  and  that  of  the  wall  to 
the  air,  respectively,  and 


(3) 


Provided  that  the  mechanical  peeling  equations  and  those  of  the 
particle  cloud  erosion  are  given,  their  effect  on'the  heat  transfer, 
shape  changes  and  thermal  stresses,  can  be  evaluated  using  the 
synthetic  analysis  computer  code. 


During  the  nosetip  ablation,  the  equation  describing  the  instant¬ 
aneous  contour  changes  is  the  following: 

where  m, ,  can  be  a  summation  of  the  mass  losses  caused  by  the  abla- 
tion,  peeling  and  erosion.  This  is  a  nonlinear  hyperbolic  type 
partial  differential  equation.  The  strategy  to  solve  this  equation 
is:  first  linearize  the  equation,  then  generate  a  difference  equa¬ 
tion  which  employs  the  implicit  difference  format  and  introduces 
artificial  viscous  terms  and,  finally,  solve  the  difference  equa¬ 
tion  using  the  iteration  method.  Details  about  these  procedures 
can  be  found  in  [5]  and  [7] .  Figure  4  which  will  be  shown  shortly 
gives  the  history  of  the  contour  changes  of  the  heat  protection 
nosetips  evaluated  using  this  algorithm  at  the  same  time  when  the 
temperature  fields  were  calculated. 


3.  INSTANTANEOUS  TEMPERATURE  FIELDS 


Since  the  carbon-base  materials  most  likely  exhibit  an  aniso¬ 
tropic  behavior,  it  is  necessary  to  calculate  the  heat  transfer 
using  the  heat  transfer  equations  for  anisotropic  materials  as  the 
starting  point.  When  computing  the  instantaneous  temperature  fields 
for  the  plug-type  nosetip  structure,  it  seems  better  to  divide  the 
domain  of  interest  into  two  parts,  the  main  part  and  the  plug  handle 
of  the  nosetipe,  between  which  an  arbitrary  boundary  such  as  a 
spherical  surface  may  be  imagined.  These  two  parts  can  be  connected 


artificially  together  as  shown  in  Figure  2.  For  the  main  part, 
after  a  spheroidal  coordinate  system  with  corresponding  origin  of 
the  coordinates  is  chosen  and  a  moving  coordinate  system  which  is 
coordinated  with  the  moving  boundary  is  introduced,  transformation 
between  the  two  coordinate  systems  is  given  as 


l  =  - 


r  —  Ti 


(5: 


Consider  the  transformation  of  the  moving  burning  boundary 
into  a  formal  fixed  boundary  and  the  heat  transfer  equations 
involving  three  principal  directions  which  are  perpendicular  to 
each  other 


ax' 

\  dx  . 

/  ay' 

<  ay  J  a::y  '  a^  / 

=  0 


:6 ) 
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Figure  2.  The  coordinate  system 
1 — artificial  boundary 

where  the  dA/d"  is  the  rate  of  thickness,  i.e.,  the  rate  of  contour 
change s, 3 in  equation  (4).  The  conditions  for  heat  transfer 
along  the  burning  boundary  are 

4  =  1.  f}t~-Q  (8) 

where 

0  =  <7.,.  +  p((^''Ccfr)(l  +  CJ  +  C;)-i|^  (9) 


Only  bodies  with  rotational  symmetry  at  zero  angle  of  attack 
and  having  isotropic  behavior  (at  least  within  the  plane  perpen¬ 
dicular  to  the  wind  axis,  i.e.,  k  =  k.  )  have  derivative  terms 

3  JD 

involving  <})  cancel.  Then  the  problem  is  simplified  to  an  axisymm- 
etric  one.  For  the  transformed  equation  (7)  as  for  an 
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axisyrranetric  problem,  the  first  and  second  derivatives  of  temper¬ 
ature  with  respect  to  the  space  coordinate  and  the  derivatives  of 
heat  flux  with  respect  to  space  coordinates  within  the  domain  all 
assume  a  central  difference  form  with  the  second  order  of  accuracy. 
Along  the  boundary,  either  the  central  difference  or  the  forward 
or  backward  difference  with  the  second  order  of  accuracy  is  assumed. 
The  derivatives  of  temperature  with  respect  to  time  take  on  the 
DuFort-Frankel  difference  form.  Then  in  the  difference  equation 
within  the  domain,  it  is  assumed  that 

T‘.  ,=  |(r;.v  +7’--:*)+o(A/‘)  ' 

where  the  superscripts  such  as  n,  denote  the  numbers  of  time  steps. 
Neglecting  the  higher  order  terms  and  rearranging  this  equation, 
one  obtains  a  three-level  explicit  difference  equation,  which  is 
suitable  for  computation  with  good  stability.  Reference  [5]  gives 
the  details  about  the  derivation  of  this  equation.  An  example  of 
an  exact  analytic  solution  using  this  algorithm  is  also  given  with 
satisfactory  accuracy  and  there  is  not  much  influence  of  the  time 
steps  on  the  computation. 


Since  there  is  no  moving  boundary  for  the  plug-handle  part, 
it  is  better  to  use  a  fixed  cylindrical  coordinate  system.  For 
rotational  bodies  at  zero  angle  of  attack  and  k^  =  kj^,  the  heat 
transfer  equation  can  be  written  in  the  following  axisymmetric 

^  f-l 

^ «  cT ,  \  3r f  /  dz  \  dz  }  (11) 


when  connecting  the  major  part  and  the  plug-handle,  one  must  con¬ 
sider  the  conditions  of  heat  flux  continuity  along  the  artificial 
boundary  between  them. 


Attention  should  be  paid  to  the  following  facts  when  solving 
equation  (11)  numerically  using  the  difference  method:  for  the  9- 
coordinate  in  the  spheroidal  coordinate  system,  the  major  part  is 
divided  into  &  uniform  grid  network.  However,  the  distribution  of 
the  grid  points  located  on  the  boundary  is  not  uniform  corresponding 
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to  direction  in  the  cylindrical  coordinate  system.  To  establish 
the  connection  conditions  of  the  normal  heat  flux  continuity  at  both 
sides  of  the  boundary,  it  is  required  that  the  grid  points  falling 
onto  the  boundary  surface  by  the  grid  network  within  the  plug- 
handle  part  coincide  with  those  of  the  main  part.  This  makes  it 
necessary  to  divide  the  plug-handle  part  into  a  non-uniform  grid 
network  to  meet  this  requirement.  In  addition,  the  boundary  is 
also  irregular.  To  harmonize  such  a  distribution  with  that  of  the 
main  part,  it  is  necessary  to  develop  a  three-level  explicit 
difference  algorithm  for  this  non-uniform  grid  network.  Specifically, 
for  the  in-domain  points,  Taylor  series  expansion  can  be  applied  to 
the  domain  using  the  temperatures  or  heat  fluxes  at  their  neighbor¬ 
ing  points  to  obtain  the  difference  approximation  to  the  temperature 
or  heat  flux  derivatives  at  those  points  with  unequal  steps.  For 
example : 
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with  the  same  technique  used  in  the  main  part,  the  three-level 
explicit  difference  equations  for  the  non-uniform  mesh  can  be 
developed . 


The  boundary  condition  for  the  irregular  boundary  is 


kt  ~  rir.  +  k,  = 

12 


(15) 
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Since  the  grid  lines  are  not  necessarily  parallel  to  the 
boundary  lines,  the  structure  of  the  difference  format  is  more 
complicated.  For  the  boundary  mesh  shown  in  Figure  3,  in  order  to 
derive  the  unequal  step  difference  expressions  to  the  first  order 
derivatives  of  the  temperatures  at  boundary  nodal  point,  (J,  1) 
with  the  second  order  of  accuracy,  the  difference  expressions  for 
the  second  order  derivatives  which  include  the  cross  terms  must  be 
derived  simultaneously  so  as  to  generate  the  difference  equations 
for  boundary  condition  equation  (15).  When  the  partial  derivatives 
are  approximated  by  the  non-uniform  difference,  there  is  no  need 
to  make  any  artificial  assumptions  usually  made  at  the  boundary, 
except  that  the  difference  is  used  to  approximate  the  derivatives. 
As  a  result,  a  more  satisfactory  solution  can  be  obtained. 

What  has  been  mentioned  above  has  been  realized  in  our  synthe¬ 
tic  analysis  program.  Calculations  show  that,  similar  to  the 
uniform  difference  format,  the  time  step  in  the  non-uniform  case 
does  not  have  much  influence  on  the  accuracy  and  stability  of 
solutions  under  the  sample  calculation  conditions. 

Results  of  the  sample  calculations  for  the  history  of  the 
graphite  plug-type  nosetip  ablation  and  heat  transfer  during 
re-entry  are  shown  in  Figure  4. 
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Figure  3.  The  mesh  around  nodal  point  (J,  1). 
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Figure  4.  Temperature  fields  and  ablation  contour  of 
the  plug-type  carbon-base  nosetip  at  different  alti¬ 
tudes  during  re-entry. 

1 — km. 


4.  STRESS  AND  STRAIN  FIELDS 


Since  the  carbon-base  materials  have  a  very  small  range  of 
plasticity,  the  objects  of  this  study  are  limited  to  those  of 
linear  elasticity  with  small  strains.  To  calculate  thermal  stress 
fields  and  thermal  strain  fields  from  the  known  temperature  dis¬ 
tribution,  the  incomplete  generalized  potential  energy  principles 
[8]  are  used,  i.e.,  the  small  displacement  linear  elasticity  prin¬ 
ciples  which  take  into  account  the  boundary  load  and  the  boundary 
displacement  conditions  and  only  satisfy  the  strain-displacement 
relationship, 

_  1 1(«<  — S#)a„n,t/s 


(16) 


with  the  finite  element  method.  Note  that  the  subscripts  i  and 
j  in  the  equation  denote  the  corresponding  components  and  the 
summation  operation.  The  axisymmetric  problems  can  be  simplified 
into  those  within  the  meridian  plane  in  the  cylindrical  coordinate 
system.  Triangles  are  selected  as  the  finite  elements,  and  in 
order  to  avoid  the  curious  integral  along  the  symmetric  axis,  refer¬ 
ence  [9]  suggests  that,  except  that  the  parameters  to  be  determined 
at  the  nodal  points  include  the  axial  displacement,  v,  the  circum¬ 
ference  strain  is  directly  assumed  as 


(17) 


A  linear  interpolating  function  can  be  constructed  in  an  element 
by  these  unknown  parameters. 

\lrUN.  IN,  m^Udy  (18) 

where  superscript,  e,  indicates  the  value  in  an  element  and  {5}® 
is  the  unknown  nodal  point  parameter  column  matrix  in  the  element. 

=Luf,  v,  w,  vj  (19) 


Use  interpolating  function  (18)  to  calculate  the  integral  in  equa¬ 
tion  (16) ,  find  the  summation  for  all  elements  to  compose  the 
general  functional.  Then  get  the  stationary  values , i . e. ,  a  set  of 
linear  equations  which  contain  a  symmetric  coefficient  matrix  in 
which  V  and  w  are  unknowns.  Then  find  the  nodal  point  parameters 
to  be  determined.  The  order  of  this  matrix  is  twice  the  total 
amount  of  nodal  points.  Accordingly,  the  strain  and  stresses  can 
be  solved  based  upon  the  relationships  between  strain  and  displace¬ 
ment  and  stress  and  strain. 


In  equation  (16),  the  boundary  conditions  are  entirely  brought 
into  the  functional.  This  not  only  makes  the  scheme  dealing  with 
the  boundary  conditions  more  reasonable,  but  also  helps  investi¬ 
gate  the  influence  of  various  types  of  boundary  conditions.  An 
urgent  problem  which  needs  to  be  solved  is  that  only  the  outside 
boundary  which  is  subjected  to  aerodynamic  pressures  has  been 
dealt  with  so  far,  while  the  more  important  contact  surfaces  have 
not  been  studied  yet,  or  only  treated  as  a  free  boundary.  In 
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order  to  fully  investigate  such  boundary  conditions,  the  author 
divides  the  boundary  conditions  into  three  categories,  then  gives 
the  derivation  and  sample  computation.  The  first  category  is  of 
free  boundary;  in  the  second  category,  displacement  on  the  boundary 
is  assigned  in  advance,  e.g.,  a  zero  boundary  displacement;  and  in 
the  third  one,  the  boundary  displacement  is  partially  assigned, 
such  as  the  case  where  only  the  tangential  displacement  is  allowed, 
but  not  allowed  in  the  normal  direction.  It  undoubtedly  deals  with 
the  mixed  boundary  conditions  involving  simultaneously  the  three 
types  of  boundaries  with  more  ease. 

5.  NUMERICAL  METHOD  FOR  SHAPE-SELECTING 

It  seems  very  clumsy  to  generate  in  advance  each  finite 
element  grid  when  computing  the  thermal  stresses  for  the  nosetips 
with  different  original  contours.  Even  for  the  nosetips  with  the 
same  original  contour,  the  contour  changes  with  time  due  to  the 
ablation.  It  is,  therefore,  almost  impossible  to  generate  each 
grid  network  manually  for  the  nosetips  with  their  instantaneous 
contours.  Thus,  thermal  stresses  can  usually  be  estimated  only 
from  a  fixed  contour  or  a  few  contour  changes.  Moreover,  each  time 
one  calculates  a  shape,  one  has  to  generate  one  grid  system.  So 
it  is  necessary  to  develop  a  special  technique  which  can  be  used 
to  generate  automatically  the  finite  element  grid  network  as  the 
shape  of  the  body  changes.  The  author  developed  such  an  automatic 
grid  generation  technique  in  editing  this  synthetic  analysis  com¬ 
puter  code.  The  finite  element  grid  mesh  generated  by  this  tech¬ 
nique  not  only  can  automatically  match  the  changing  shapes  caused 
by  the  ablation,  but  also  can  automatically  match  various  original 
contours  with  different  dimensions,  provided  that  the  adjustable 
dimensions  of  the  original  shapes  are  available.  As  shown  in  Figure 
5,  after  the  adjustable  parameters  R  ,  z  ,  z  ,  3  and  0^  are  pro- 
perly  selected,  the  finite  element  can  be  computed  for  various 
bodies  with  the  required  shapes.  The  shape-selection  computation 
can  be  easily  performed  after  selecting  the  different  combinations 
of  these  adjustable  parameters.  Furthermore,  combining  the 


Fig.  5.  The  adjustable  paran.eters  of  the  initial  contour  of  the 
pi  Uf,- type  nose  tips 


Fig.  C.  Flow  chart  of  tile  compute  code 

1-  Criminal  contour  and  temperature  field  at  the  re¬ 
entry  peint;  ?-  Trajectory  time,  t;  3-  Difference  mesh; 

TicT-rr.ai  envii-(;nir.enl  ;  h-  Instantaneous  temperature 
field;  t~  Ablation;  7-  boundary  conditions;  B-  Rigidity 
t.atrix  ind  load  vector;  ‘■j-  Finitt.'  element  mesh; 


IF-  1 nc tan taneous  contour;  11-  Stress  and  strain; 
1j-  si  sp  j  ac;ef.'ien  L  . 


quantitative  analysis  of  the  effects  of  different  boundary  condi¬ 
tions,  it  is  possible  to  provide  a  useful  basis  for  the  engineering 
designs . 

The  synthetic  analysis  computer  code  edited  starts  the  cal¬ 
culations  from  the  re-entry  point,  then  for  the  thermal  environ¬ 
ment,  based  on  the  trajectory  and  body  parameters,  followed  by  the 
estimation  of  the  ablation  and  instantaneous  temperature  fields. 

The  calculations  for  the  thermal  environment,  ablation  and  instant¬ 
aneous  temperature  fields  are  coupled.  The  code  computes  the 
instantaneous  contours,  based  on  the  ablation  calculation  results; 
from  the  instantaneous  contour  and  instantaneous  temperature  fields 
and  considering  various  types  of  boundary  conditions,  the  code 
calculates  the  stress  fields  and  the  strain  fields  at  that  instant. 
The  difference  grid  mesh  and  the  finite  element  grid  mesh,  changing 
with  the  instantaneous  contours,  are  automatically  generated  at 
every  instant.  A  brief  flow  chart  of  the  complete  computer  code 
is  shown  in  Figure  6. 

Since  the  time  step  can  be  taken  as  long  as  0.2  second  when 
computing  the  temperatures  and  the  contours,  the  time  spent  during 
computing  all  these  parameters  along  a  trajectory  is  limited. 

As  an  example,  results  of  the  selecting  calculations  of  the 
graphite  nosetips  for  one  of  the  parameters  are  shown  in  Figure  7 

6.  CONCLUSIONS 

A  synthetic  analysis  algorithm  in  which  the  thermal  environ¬ 
ment,  ablation,  heat  transfer  and  thermal  stresses  are  considered 
and  an  automatic  grid  generation  technique  has  been  developed,  is 
briefly  introduced.  With  this  synthetic  analysis  computer  code, 
it  is  possible  to  perform  a  comprehensive  numerical  analysis  for 
ablation,  heat  transfer  and  thermal  stresses  during  re-entry  of 
the  carbon-base  nosetips  into  the  atmosphere.  Calculations  were 
made  for  nosetip  shapes  and  the  connection  format  of  their  contact 
surfaces. 
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Abstract 

In  this  paper,  the  principle  of  synthetic  analysis  calculation  and  some  treatment  technics 
about  reentry  thermal  environment,  ablation,  conduction  and  thermal  stresses  of  vehic¬ 
les  are  described  briefly.  The  thermochemical  ablation  calculation  for  chemical  reactions  on 
the  ablating  surface  of  carbonbase  materials.  is  simplified  by  using  selective  calculation 
result.  The  temperature  fields  are  calculated  by  means  of  the  three  level  explicit  difference 
schemes  for  both  constant  and  unconstant  space  steps  within  total  region.  The  thermos- 
tresses  are  calculated  by  the  finite  element  method,  which  is  based  on  the  noncomplete 
generalized  potential  energy  variational  principles  under  the  supposed  condition  of  linear 
elasticity.  The  calculation  method  of  the  effects  for  various  displacement  boundary  con¬ 
dition  are  provided.  An  automatic  mesh  dividing  technic  is  realized  in  the  comprehensive 
analysis  procedure  and  selective  calculations  for  shapes  can  be  made  with  this  technic 
conveniently. 
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